function [NetEst_Output] = Fn_EstimateNetwork_Mar2022(starts, Cons, G_ij, est_tol, DELTA)
% Inputs: Covariates (Cons), Nets and Outcomes (NetOut)
% Settings: starts, est_tol
% Output: 
% Note that this has option DELTA for delta-missingness method. Set DELTA =
% 0 for normal method

% Set bounds for estimation
Lb = -inf*ones(4*Cons.dimX+3+Cons.n,1);
Lb(1) = 0;
Lb(2+Cons.dimX) = 0;

Ub = inf*ones(4*Cons.dimX+3+Cons.n,1);
Ub(1) = 1;


%% Estimate
% Replace G_ij when missing
G_ij_est = G_ij - Cons.NETMISS*DELTA;

options = optimoptions('fmincon','DerivativeCheck','off','Display','iter-detailed','GradObj','on','TolFun',est_tol,'maxiter',10000,'MaxFunEvals',10000);
fun1 = @(x)GMM_Fn_Mar2022(x,Cons,G_ij_est,[]);
[BGDA_hat, fval_min, exitflag] = fmincon(fun1,starts,[],[],[],[],Lb,Ub,[],options);


BGD_hat = BGDA_hat(1:4*Cons.dimX+3);
A_hat = BGDA_hat(4*Cons.dimX+3+1:4*Cons.dimX+3+Cons.n);


%% To get V_C
A_i = Cons.I_i*A_hat;
A_j = Cons.I_j*A_hat;
%A_j_doti = Cons.Trans_doti*A_j;

G_ji = Cons.Flip_ij*G_ij;
%G_ij_doti = Cons.Trans_doti*G_ij;
%G_ji_doti = Cons.Trans_doti*G_ji;

RHS = [G_ji, Cons.X_j, A_j, bsxfun(@times,Cons.X_i,Cons.X_j), bsxfun(@times,Cons.X_i,A_j), bsxfun(@times,A_i,Cons.X_j), bsxfun(@times,A_i,A_j)];
C_hat_raw = -G_ij + RHS*BGD_hat + Cons.NETMISS*DELTA;
lambda_hat = inv(Cons.I_i'*Cons.I_i)*Cons.I_i'*C_hat_raw;
C_ij_hat = C_hat_raw - Cons.I_i*lambda_hat;

if sum(Cons.NETMISS) < size(Cons.NETMISS,1)
    V_C = bsxfun(@times,ones(size(Cons.NETMISS,1),1) - Cons.NETMISS,[C_ij_hat, Cons.Flip_ij*C_ij_hat])'*bsxfun(@times,ones(size(Cons.NETMISS,1),1) - Cons.NETMISS,[C_ij_hat, Cons.Flip_ij*C_ij_hat])/(size(Cons.NETMISS,1)-sum(Cons.NETMISS));
else
    V_C = [C_ij_hat, Cons.Flip_ij*C_ij_hat]'*[C_ij_hat, Cons.Flip_ij*C_ij_hat]/size(Cons.NETMISS,1);
end


%% To get M_hat and beta_M_hat
l_M_hat = log(Cons.I_i'*(exp(G_ij).*exp(C_ij_hat)));

X_M = [ones(size(Cons.X_raw,1),1), Cons.X_raw, Cons.P, bsxfun(@times,Cons.X_raw,Cons.P)];

beta_M_hat = inv(X_M'*X_M)*X_M'*l_M_hat;

sigma2_M_hat = mean((l_M_hat - X_M*beta_M_hat).^2);

NetEst_Output = struct('BGD_hat',BGD_hat,'A_hat',A_hat,'V_C',V_C,'l_M_hat',l_M_hat,'beta_M_hat',beta_M_hat,'sigma2_M_hat',sigma2_M_hat,'fval_min',fval_min,'exitflag',exitflag);

end

